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SUMMARY 


This paper considers complex transcendental eigenvalue problems where one is interested in pairs 
of eigenvalues that arc restricted to take real values only. Such eigenvalue problems arise in dy- 
namic stability analysis of non-conservative physical systems - flutter analysis of aeroelastic systems, 
to name one example. Some available solution methods are discussed and a new r method is pre- 
sented. Two computational approaches are described for analytical evaluation of the sensitivities 
of these eigenvalues when they arc dependent on other parameters. The algorithms presented are 
illustrated through examples. 


♦NASA Resident Research Associate. 


Introduction 


Dynamic analysis of physical systems often requires the solution of an eigenvalue problem. Opti- 
mization of dynamic systems also requires efficient computation of the derivatives of the eigenvalues 
and eigenvectors. The computation of these derivatives is referred to as sensitivity analysis. De- 
velopment of numerical algorithms for the evaluation of eigenvalues and eigenvectors (Wilkinson, 
1965) and, recently, their derivatives with respect to system parameters (e.g., Murthy and Haftka, 
1988), for the linear algebraic eigenvalue problem 

(A-;j)u = o (i) 

where the matrix A is complex non-hermitian has attracted extensive research efforts. The 
eigenvalues X and eigenvectors u for such problems are also generally complex. Some of these 
methods and algorithms can be directly extended to the ^-matrix (Frazer, Duncan and Collar, 1960) 
eigenproblem 

A(;.)ti = 0 (2) 

where the elements of the complex non-hermitian matrix A are polynomials in the eigenvalue, /. 
A discussion of the /-matrix eigenproblem can be found in Rokne (1985). Extensive analytical 
results on the sensitivities of the eigenvalues and eigenvectors of eq. (2) are presented by Taylor and 
Kane(1975) for the special case where the polynomials are quadratic. These sensitivities were used 
by Fritzen and Nordmann (1982) for predicting the stability behavior of non-conservative rotor 
dynamic models. 

In this paper, we are concerned with the solution and sensitivity analysis of the eigenvalue 
problem 


a(/i,;. 2 ) u = o (3) 

where the elements of the complex non-hermitian matrix A are transcendental functions of the pair, 
A, and A 2 , which are restricted to take real values only. Comparing eq. (3) to the A-matrix 
eigenproblem of eq. (2), we note two differences: 1) the pair of real numbers A, and A 2 replace the 



complex number / and hence play the role of eigenvalues and 2) the elements of A in eq. (3) are 
transcendental functions of the eigenvalues instead of polynomial functions. We also note that, if 
and A 2 are allowed to take on arbitrary complex values, eq. (3) does not represent a well-defined 
eigenvalue problem. Hie eigenvectors u in eq. (3) remain complex-valued. 

This problem was motivated by the need to automate the flutter analysis in acroclastic opti- 
mization problems. Murthy and Kaza (1988) automated the flutter analysis procedure using the 
direct solution approach which results in the eigenvalue problem of the form given by eq. (3). The 
flutter analysis procedures of Cardani and \lantegazza(1978) and Meyer(1988) also result in 
eigenvalue problems of this form. Here, the matrix A involves the generalized stiffness, mass and 
aerodynamic matrices of the aeroelastic structure, which is undergoing steady-state oscillations in 
a state of neutral stability. The generalized stiffness, mass and aerodynamic matrices are usually 
computed by specialized structural and aerodynamic analysis programs which are computationally 
intensive. Depending upon the formulation, the eigenvalue Ay is a speed parameter such as the 
Mach number or the rotational speed at flutter and the eigenvalue A 2 is a frequency parameter such 
as the reduced frequency at flutter. In aeroelastic literature, the pair A ] and A 2 is called the matched 
flutter point and the eigenvector u the flutter mode. The viewing of the flutter analysis problem 
as a complex eigenvalue problem involving a pair of real eigenvalues was first proposed by Frazer 
(1946) and is called the direct solution approach. An alternate approach to flutter analysis, known 
as the U-g method (Bisplinghoff, Ashley and Halfman, 1955), which decomposes the solution of 
eq. (3) into a series of linear algebraic eigenvalue problems, has been more popular in aeroelastic 
literature, perhaps because of the increased physical insight it offers. However, the direct solution 
approach is more attractive in terms of automating the flutter analysis procedure and also more 
efficient, as discussed by Murthy and Kaza(1988), because it avoids troublesome eigenvalue tracking 
and replaces the inner-outer iteration loop required by the U-g method by a single iteration loop. 

The objective of the present paper is to present (i) a solution method for calculating the 
eigenvalues and the eigenvectors of eq. (3) and (ii) methods for sensitivity analysis of eq. (3). The 
solution method was developed by Murthy and Kaza (1988) but is presented here in a more general 
context and with additional examples. The sensitivity analysis methods appear to be new. It is 
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assumed that the computation of the matrix A is much more expensive in terms of CPU time than 
its inversion. This is typically the case in aeroelastic applications (Murthy and Kaza, 1988). 


Solution Methods 

The eigenvalue problem of eq. (3) differs from those of eqs. (1) and (2) in many respects. While 
the conventional eigenproblems of eqs. (1) and (2) define eigenpairs, the eigenvalue problem of eq. 
(3) defines eigentriples, each consisting of A j , A 2 and u. Also, it is well-known that the number of 
the eigenpairs of eqs. (1) depends only on the order of the matrix and, in the case of eq. (2), on the 
matrix order and the highest degree of the polynomials in A. In the case of eq. (3), the number of 
the eigentriples is not known in advance and could be infinite, finite or none. In addition, while 
the eigenvalue problems of eqs. (1) and (2) admit bi-orthogonal relationships between the left and 
the right eigenvectors and can be reduced to standard canonical forms, no such relationships and 
canonical forms are known to exist for the eigenvalue problem defined by eq. (3). This has im- 
portant consequences in terms of sensitivity analysis and is discussed further in the next section. 

In this paper, we will restrict our attention to those problems that satisfy the following as- 
sumptions: 

• There exist real values of k } and A 2 such that eq. (3) is satisfied for some vector u. 

• For such values of and A 2 , there exists a neighborhood in the A 1? A 2 -plane in which the 
members of the matrix A are continuously differentiable functions of Aj and A 2 . 

These assumptions are usually satisfied in aeroelastic applications. Solution of eq. (3) consists 
of finding the pairs of values of A } and A 2 such that eq. (3) is satisfied for a non-zero vector u, and 
the determination of the vector u itself within a constant complex multiplier. 

There are essentially two different approaches to the development of numerical schemes for 
this problem. In the first approach, the condition of a vanishing determinant for the existence of 
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a non-zero vector u is used to obtain the values of A, and A 2 . The eigenvector u is obtained inde- 
pendently afterwards. Thus, 

det[A(A,, A 2 )] = A 2 ) = 0 (4) 

Eq. (4) is one complex equation in two real unknowns, A, and A 2 , and can be solved by any of a 
variety of quasi-Newton methods. Eq. (4) is first rewritten in terms of its real and imaginary parts 
as 


0r(A i, A 2 ) = 0 
^/(Aj, A 2 ) = 0 


Eqs. (5) represent two nonlinear equations in two unknowns, A, and A 2 . The typical iteration 
scheme for their solution would then be 





( 6 ) 


where k is the iteration number and J is the Jacobian matrix of eq. (5) and is given by 

5 _£_R 5 _£r_ 

d/.j <?A 2 

J = dDj dDj 

H7 ~dJJ 


(?) 


If the Jacobian is analytically evaluated, eq. (6) represents Newton's method. Quasi-Newton 
methods approximate the Jacobian in some manner. Hassig (1971) and Stark (1984), for example, 
have used this approach, solving eq. (4) by the Regula Falsi and the secant methods respectively. 


In the second approach, proposed by Cardani and Mantegazza (1978), the eigenvalue prob- 
lem of eq. (3) along with a normalizing condition on the eigenvector u are treated as a set of non- 
linear equations in the unknowns, u, /j and / 2 . These nonlinear equations are then solved using 
standard techniques, obtaining the complete eigentriple simultaneously. While this approach ob- 
viously results in a much larger system of non-linear equations than the previous determinant 
equation approach, it has recently been advocated by Mever (1988) over the determinant equation 
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approach because “the basic drawback to solving the determinant equation is that since it is im- 
practical to compute derivatives of the determinant with respect to system parameters, only first- 
order convergent methods such as Regula Falsi or a secant method can be used, instead of more 
rapidly-converging methods such as Newton's method ” and “ convergence to the wrong root, 
known as mode switching, can occur when two modes are nearly equal Meyer (1988) presented 
a continuation method for the solution of eq. (3) using the second approach. 

This paper addresses the first drawback noted by Meyer (1988). A practical iterative scheme, 
for solving the determinant equation with a convergence rate that is more rapid than the secant 
method and approaching that of Newton's method, is proposed. It is first shown that all the de- 
rivatives of the determinant can be evaluated at the cost of a single matrix inversion. Further, a 
quasi-Newton method, based on Broyden's updates for the derivatives of the matrix A (rather than 
for the Jacobian matrix as is done in Broyden's method), is shown to be more rapidly convergent 
by means of examples. The second drawback of converging to the wrong root will not be addressed 
in this paper. However, this is felt to be only a minor limitation because most practical aeroelastic 
structures do not often give rise to nearly equal modes that flutter. 

The eigenvector u (flutter mode, in aeroelastic literature) and the left eigenvector v can be 
computed very efficiently using the information available at the time of convergence of the iterative 
scheme in an inverse iteration procedure. 


Computation of Eigenvalues 

We consider the computation of the eigenvalue-pairs for eq. (3) by solving eqs. (5) through the it- 
erative scheme of eq. (6). Note that analytical evaluation of the Jacobian requires the derivatives 
of the characteristic determinant, D . The derivatives of the characteristic determinant are computed 
by using the following result, known as the Trace theorem, proved by Lancaster(1964): 

If the elements of A(a) are differentiable functions of a, then for any * for which the charac- 
teristic determinant £>(a) ^ 0, 

= D(=<) . trace of (A -1 ) (8) 
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Using eqs. (7) and (8), the iterative scheme of eq. (6) can be written as 


M _JM i m 

A 2 j^ +1 ) {;. 2 J W (gRhj - h R gj\ k) \-g] j<*) 


(9) 


where 


gR + *gl = trace of 
Hr + ih] = trace of 




( 10 ) 


Thus, the trace theorem gives us a practical method for calculating the derivatives of the de- 
terminant. Note that, in evaluating eq. (10), matrix multiplication need not be performed, since 
only the trace is required. The computational cost of the matrix inversion required in eqs. (10) is 
not significant if the use of such derivatives improves the convergence rate of the determinant 
equation solution. This is because the cost of inverting the matrix A is usually negligible in com- 
parison to the cost of evaluating it. The iteration is terminated when Aj and / 2 satisfy some con- 
vergence criteria or the matrix A cannot be inverted. 

In eq. (10), if the derivatives of the matrix A with respect to the eigenvalues A l and ). 2 are 
evaluated analytically, then eq. (9) is equivalent to Newton's method and thus gives quadratic 
convergence. Analytical differentiation of A is not often practical and one resorts to approximate 
evaluation of the derivatives of A, sacrificing quadratic convergence and hence increasing the num- 
ber of iterations required for convergence. Approximate derivatives of A can be evaluated simply 
by finite-differences. However, finite-differencing increases the computational cost of each iteration, 
further degrading the efficiency of the iterative procedure. It is therefore proposed that the deriva- 
tives of A be evaluated during the iteration using Broyden's update formula (Dennis and More, 
1977), which is used in Broyden's method to update the Jacobian. Numerical examples shown later 
indicate that this updating scheme, in combination with eqs. (9) and (10), results in an iterative 
scheme with a high rate of convergence (close to that of Newton's method) with the same cose per 
iteration as that of secant and Broyden's methods. 
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The elements of A are assumed to be linear functions of A, and A 2 in the direction of the last 
iterative move. Hence, the derivatives of A at the A:-th iterative solution satisfy 

A (*-1) = A (*) + ('f^)(*) (A l(*-l)”‘ i l(*)> + (i^)w( A 2(*-l)-^)) (11) 


In addition, the derivatives of A in the direction orthogonal to the last iterative move are as- 
sumed to be the same at the A-th and the (k — l)-th iterative solutions. This implies 


( "Ja^ ~ + ( )w (A 2 ~ A 2(k)) 


= (~d^ + ( )(k-l)( Jl 2 - ; -2(A')) (12) 


for any Aj and A 2 satisfying 


( ; - 1 ~ 1 (A)) ( ; - 1 (A- 1 ) " ; T (*)) + ( ; 2 ~ h(k)> ( ; -2 (k- 1 ) “ h (A)) = 0 


(13) 


Eqs. (11)-(13) uniquely determine ^~^~^(k) • Esing the 

6' l -\(k) = ^-l(A-l) ~ f - 1 (k) ^ b)-2(k) - *2(k-\) - *2(k)< ^(*) * s 


notation 


obtained as 


(^-)w 


< A (*-1)- A (*))' S; 1(«:)+ (^)(i-l) W 2(*)- (^)(*-l)^]Wj ,S; '2(<r) 


. 2 2 

6> \(k) + 6 h{k) 


•( 14 ) 


A similar expression can be obtained for \k) by simply exchanging the subscripts 1 

and 2 for A in eq. (14). 


Computation of Eigenvectors 
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The eigenvector u (flutter mode, in aeroelastic literature) can be computed very efficiently using the 
information available at the time of convergence of the iterative scheme, eq. (9), in an inverse iter- 


ation procedure. If u (0) is any vector that is not totally deficient in the eigenvector, the following 
inverse iteration scheme can be used to calculate the eigenvector. 

u (k ) = A lu (A-l) k= 1,2,... (15) 

At the time of convergence of the iteration given by eq. (9), the matrix is already available in 
decomposed form and can be used directly in eq. (15). 

Sensitivity calculations sometimes also need the left eigenvector associated with and ). 2 . 
The left eigenvector is defined by the equation 

v 7 a(;. 1 ,;. 2 ) = o 7 ' (16) 

Once again, inverse iteration can be used to calculate the left eigenvector as follows. If v^ 0) is any 
vector that is not totally deficient in the left eigenvector, 

y (k) ~ [ A? 3 S*— 1) k= 1,2,... (17) 

As in the case of the linear eigenvalue problem, one step of the inverse iteration scheme is normally 
sufficient to obtain the eigenvector within a reasonable accuracy. There is again no need for another 
decomposition to solve eq. (17), as [A 7 ^] -1 = [A^J^and the previous decomposition can be used 
by replacing backward substitution with forward substitution. 


Numerical Examples 

The performance of the above algorithms is demonstrated through two numerical examples. Both 
examples present complex transcendental eigenvalue problems having known pairs of real 
eigenvalues. Example 1 uses trigonometric functions and example 2 uses hyperbolic functions. 
The examples were designed to illustrate the various features of the algorithms and for easy verifi- 
cation of the results. 
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Example l:For this example, the complex transcendental matrix is given by 


... . . pi MPih + Pl) 

-1 ’ 2 z 2 sin(ajAi + a 2 ) 


( 18 ) 


This eigenproblem has an infinite number of pairs of real eigenvalues. A family of such pairs is 
given by 


;-i = 


m — <t 2 


sn-P 2 
2 ~ h 


r= 0, ±1,±2, ... 5 = 0,± 1,±2,... 


(19) 


The nominal design point is represented by the values 
Z\ = 1 + i z 2 = 1 - / 

«1 = 2 = — 1 ( 20 ) 
fil=2 fi 2 = -3 

The eigenvalue pair of interest is given by 2] =0.5 and k 2 = 1.5. 


Example 2:In the second example, the eigenproblem involves the complex transcendental matrix 
A(/ l5 2 2 ), the elements of which are given by 


a pq = z pq sinh(a p ;., + P p ) + zj$ sinh(a 2 / 2 + S p ) 


( 21 ) 


p= 1,2, ... ,n q= 1, 2, ... ,n 


ztyj and zfij represent complex constants chosen randomly for this example. <x p , (i p and S p are real 
quantities that serve as design parameters. It can easily be shown that when the elements of the 
matrix A are given by eq. (21), the eigenproblem of eq. (3) has at least n pairs of real eigenvalues 
that are independent of the constants zQJ and These are given by 


;-i = - 




P = 1,2, ... , n 


( 22 ) 


10 



There are additional pairs of real eigenvalues that depend on z^ and z^. 


The results presented here for example 2 are for n = 6 and the nominal design point is re- 
presented by the values 

*P = P 

P p = —0.05 p 2 (23) 

6 = -0.005/? 3 
for p & 2 and 
“ 2 = 1 

h = -° 5 ( 24 ) 

S 2 = — 0.1 

The eigenvalue pair of interest is given by / j = 0.5 and / 2 = 0.1. 


Rate of convergence 

The rate of convergence of the present solution procedure was earlier shown by Murthy and Kaza 
(1988) to be more rapid than that of the generalized secant method for the case of the aeroelastic 
analysis of propfan blades. In this paper, further results are presented comparing the present sol- 
ution procedure, represented by the iterative scheme of eq. (9), to Broyden's method and Newton's 
method. In Broyden's method, the directional update is applied to the Jacobian matrix. The 
present procedure differs from Broyden's method in that the directional update is applied to the 
system matrix rather than the Jacobian matrix. The eigenproblems associated with both the above 
examples are solved by Broyden's method, the present procedure and Newton's method. For 
Broyden's method and the present procedure, the exact analytical derivatives were used to calculate 
the initial Jacobian matrix. For Broyden's method and the present procedure, only the matrix A 
has to be evaluated once per iteration. For Newton's method, the matrix A and its analytical de- 


ll 



rivatives with respect to Aj and A 2 are evaluated once per iteration. The same convergence criterion 
was used for all results. Tables 1 and 2 show, for the eigenproblems posed by examples 1 and 2 
respectively, the number of iterations required for convergence starting with various initial guesses. 
The tables demonstrate that the present procedure is far superior to Broyden s method in terms of 
the number of iterations required for convergence. Comparison to Newton's method also demon- 
strates that at most two more iterations are required with the present procedure than with Newton's 
method. Thus, the present solution procedure appears to have a rate of convergence that is close 
to that of the quadratically convergent Newton's method. 


Sensitivity Analysis 

Let the matrix A and hence the cigentriple A } , A 2 and u be functions of design parameter vector p 
with individual parameters denoted by Greek subscripts, e.g. p a . We assume that the design pa* 
rameters are all real and that the pair A x and ?. 2 * s a distinct solution for eq. (4). 


Eigenvalue Sensitivities 


The objective is to obtain expressions for the eigenvalue sensitivities 
tiating eq. (4) with respect to the design parameter p a , we have, 


3 A | 

'~dP* 


and 


3 A 2 


Differen- 


dP 3D dA ] dD d _h _ 0 

dp a + dAj dp a dk 2 dp a 


(25) 


Eq. (25), once again involves the derivatives of the determinant of A , this time evaluated at the 
solution obtained at the convergence of the iterative scheme given by eq. (9). However, the trace 
theorem in the form of eq. (8) cannot be used to calculate these derivatives, as the the determinant 
vanishes at the solution A } and / 2 - Hence, eq. (8) is rewritten as 


= trace of (adj\ ) 
3a y J 3a } 


(26) 


using the relation between the inverse and the adjoint matrices, 
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( 27 ) 


A-i ad i x 
A D 

Direct calculation of the adjoint matrix for use in eq. (26) is computationally expensive and can 
be avoided as shown below. 

Frazer. Duncan and Collar (1960) have shown that, if A(A) is a X -matrix and X 0 is a simple root 
of the characteristic equation det A(2) = 0, then the rows and columns of the adjoint matrix, 
adjA(X 0 ) , are directly proportional respectively to the left and right eigenvectors corresponding to 
Xq. It can similarly be show n that the adjoint matrix has the same property in the case of the 
transcendental eigenprobem with real eigenvalues, given by eq. (3), if the left eigenvector is defined 
by eq. (16). Thus, if X x and X 2 form a simple root pair of eq. (4) and u and v satisfy eq. (3) and 
eq. (16) respectively, then adjA(X x , X 2 ) has rows and columns directly proportional to v and u re- 
spectively. That is, 

adjA(Xy, X 2 ) = du\ T (28) 

where d is some non-zero constant dependent on the eigenvector normalization. Thus, the deriv- 
atives of the determinant at the solution, X ] and X 2 , can be expressed as 

dD j , f / T dA v 

— — = d trace of (uv -r— ) 

Ox Ox 


Or, 


dD . T d A 


u 


(29) 


Using eq. (29) to evaluate the derivatives of the determinant in eq. (25) and cancelling out the un- 
known constant d, we obtain the sensitivity equation 


T d A , T d.\ d/ - 

V — U + V u 


d P a 




1 °Pa 


1 + ,TJd. a ik . 0 




d Px 


(30) 


Since the eigenvalues X } and X 2 are restricted to real values and the design parameters are assumed 


a;., 


dX n 


to be real, the sensitivities of the eigenvalues, — — and 

d P* We 


must also be real. Thus, eq. (30) 
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<5/1, 


<5/U 


is a complex equation, in the two real unknowns, and . Unlike eq. (4), the sensitivity 
equation is linear and it can be solved exactly. Sobieski(1986) pointed out that the sensitivity 
equation is always linear, even for nonlinear systems. Solution of eq. (30) gives us the following 
expressions for the eigenvalue sensitivities 


a;., 


a;, 2 


J m(y T M- U vTj£_ jj) 


d P a 


, , T d A -T d A -s 

Im{y aI7 uv alj^ 


i m(y TJ± u ,rj± d 






a;., 


(31) 


(32) 


where the bar indicates the complex conjugate. These expressions are essentially the discrete 
equivalents of the flutter point sensitivities developed by Seyranian (1982) for the case of a slender 
elastic beam in incompressible flow. Note that, while the eigenvectors are only determined up to 
an arbitrary complex multiplier in the absence of appropriate normalizing conditions, the eigenvalue 
sensitivities are independent of the normalizing condition. 


Integrated Computation of Eigenvalue and Eigenvector Sensitivities 

Murthy and Haftka (1988) described two computational approaches, the Direct approach and the 
Adjoint approach, for calculating the sensitivities of the eigenvalues and eigenvectors of linear al- 
gebraic eigenvalue problems and showed that both the approaches are competitive, under different 
sets of conditions. The Adjoint approach requires all the left and the right eigenvectors and the 
Direct approach requires no left eigenvectors. The lack of bi-orthogonality relationships in the case 
of the eigenvalue problem of eq. (3) eliminates the possibility of utilizing the Adjoint approach for 
calculating derivatives of the eigenvectors, because the bi-orthogonality property is crucial for such 
an approach. 

When eigenvector sensitivities are required, both the eigenvalue and the eigenvector sensitiv- 
ities can be computed in an integrated manner, using the Direct approach. Bindolino and 
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Mantegazza (1987) presented one such method, which operates in the real domain rather than the 
complex domain. The normalization condition they used for the eigenvectors is unreliable when 
the eigenvectors arc complex. Murthy and Haftka (1988) discussed the merits of different normal- 
ization conditions to be used in sensitivity analysis computations. In the following, a Direct ap- 
proach, requiring no left eigenvectors and based on proper eigenvector normalization, is developed 
for the integrated computation of the eigenvalue and eigenvector sensitivities for the eigenproblem 
posed by eq. (3). 


Direct differentiation of eq. (3) with respect to the parameter p a gives 


5 A , d\ d h , 0A a; -2 

-t U -f “T^ — U — + -r: — U — — 

dp a a/i a p a d/. 2 dp u 


+ A 


du 

d P* 


= 0 


which can be w'ritten in matrix form as 


(33) 



(34) 


Eq. (34) consists of n linear equations in (n + 2) unknowns, w 7 hich are the 

0 /] d /-2 

— — and — — and — — . Thus, two additional conditions must be stipulated 

op* d P* °p* 

a solution of eq. (34). One of these is obtained by normalizing the eigenvector 


Sky 




and the other by restricting — — and to real values. 


d P* 


d P* 


n components of 
in order to obtain 
to make it unique 


The first condition is easily enforced by normalizing the eigenvector so that a non-zero ele- 
ment is unity. Let u m be such an element. Then, since u m — 1, 


du m 

d P* 


= 0 


(35) 


and the m - th column of the coefficient matrix in eq. (34) can be deleted. Eq. (35) also reduces the 
number of unknowns in eq. (34) by one. 
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The second condition is best enforced by partitioning eq. (34) as follows, 
equation of the system given by eq. (34) separately, we have 



+ “1 


a;., d?. 2 

d P a + " 2 d P* ~ 


n 


a 


and 


dx , 3; lj dk 


where 


C = A 


/ — th row and m— th column deleted 


x = u 


m — tli element deleted 




th element deleted 


th element deleted 




th element deleted 


a 1 = / — th row of A with m — th column deleted 


? 7 l = / — th element of 


y] 2 ~ l —Xh element of 


» 7 a = / — th element of 


{*■} 

{*■} 

It”} 


dx 

Eliminating from eqs. (36) and (37), we have 


Writing the /-th 


(36) 


(37) 



(38) 


dX , 


dX- 


fal “ b / n l) ~aT- + (>J2 ~ b / n 2) IT- = ~ (»f« 






b/"«) 


where 


b / =[c r r 1 a / 


The condition of real eigenvalue sensitivities is now enforced on eq. (38) to obtain 


5£i 


^[(>7a ~ b/wgK^2 ~ b/ng)] 
M(»/l ~ b/ n,)(^ 2 ~ b/n 2 )] 


(39) 


^2 

d P* 


ImUVx ~ b/ r n a )(> ?) - b/njj] 
M(>f2 ~ b/ 7 n 2 )(^i - b/n,)] 


(40) 


Thus, the following procedure can be used to obtain the derivatives 


<?u 

3ft, 1 


and 


5/ 2 


1. Normalize the eigenvector such that u m — 1. 


2. Form the LU decomposition of the matrix C. 


3. Solve b/= by forward substitution. 


4. Calculate ■— and — from eqs. (39) and (40). 
®Py. cp a 


5. Calculate from eq. (36) bv backward substitution. 

d P* 

s TT A d\ . c’u 1 44 . CU m 

6. Expand — — to — — by settmg — — — 0 . 

C'Pa dp* dp a 

The proper choice of the two indices / and m is of great importance for the above procedure 
to calculate accurate sensitivities. The indices must be chosen so that the resulting matrix C is 
well-conditioned. Note that C is obtained from the singular matrix A(2 ]} / 2 ) by deleting the row 
corresponding to the index / and the column corresponding to the index m. Hence, in the absence 
of multiple solutions, for the matrix C to be non-singular, the /-th row and the m-th column of A 
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must be linearly dependent on the other rows and the other columns respectively. In other words, 
C is non-singular if and only if 0 and u m ^ 0. (For a detailed discussion of this point, see 
Nelson, 1976). Therefore, to obtain a well -conditioned matrix C, a reasonable choice of the indices 
/ and m would be such that v, is the element of largest absolute value in v and u m is the element of 
largest absolute value in u. Thus, while the direct method does not require the left eigenvector in 
principle, it may still be necessary to calculate the left eigenvector in order to choose the index / that 
results in a well-conditioned matrix C . 


Finite Difference Sensitivity Analysis 

The eigenvalue and eigenvector sensitivities can also be calculated by finite differences. While the 
finite difference approach has the virtues of simplicity and ease of implementation, it is known to 
suffer from efficiency and accuracy problems. For example, the finite difference approach requires 
the solution of the nonlinear equation, eq. (3), once for each design variable. On the other hand, 
the analytical approach requires the solution of only a linear equation, eq. (30) or (33), with the 
same coefficient matrix for all design variables. However, Haftka (1985) has pointed out that, for 
iteratively solved problems, the finite difference methods are usually superior to the analytical 
methods in terms of efficiency, as the solution at the nominal point is usually a very good initial 
guess for the solution at the perturbed point. 

The accuracy problems still remain, because of the truncation and the round-off errors in- 
herent to the finite difference approach. For iteratively solved problems, these difficulties are 
exacerbated due to the presence of errors arising from the early termination of the iterative process 
and the different initial guesses used for the nominal and the perturbed solutions. Haftka (1985) 
presented a novel technique for controlling the magnitude of the errors associated with iteratively 
solved problems. More accurate derivatives are obtained by treating the converged solution as the 
exact solution of an approximate equation instead of the approximate solution of an exact equation. 
The finite difference approach is not considered in this paper but is mentioned here to alert the 
reader to an important alternative to the analytical approach. 
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Numerical Examples 


A computer program was developed that implements both of the approaches to the sensitivity 
analysis. The above algorithms were applied to the two numerical examples considered earlier. 
Both examples present complex transcendental eigenvalue problems having known real valued 
eigenvalues that are dependent on certain parameters in a known manner. 


Example 3:For this example, sensitivity analysis of the eigenproblem presented in example 1, with 
respect to the parameters a ] and , is considered. The analytical sensitivities of the family of 
eigenvalue-pairs given by eq. (19) can be expressed as 


rn — <*2 ^2 q 

a “l ~~ <xj da l 


(41) 


dfii 


= 0 


d/i sn — /?2 

^7 = 


(42) 


r = 0, ± 1, ± 2, ... 5 = 0, ± 1,±2, ... 


Example 4:For this example, the eigenproblem is the one described in example 2 and the design 
parameters considered are a„, fi p and S p . The sensitivities of the eigenvalues given by eq. (22) can 
analytically be expressed as 


Pp d?.2 ^< 5 p 



(43) 


d'-i 1 

dp p a P 



(44) 



= 0 


a;. 2 i 



(45) 


The algorithms presented above correctly calculated the sensitivities of the pairs of real eigenvalues 
for the example problems. 
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Inaccuracy of the derivatives calculated by eq. (14) 


Eq. (31) and (32) require the derivatives of the matrix A with respect to the design parameter of 
interest, p a , and the two eigenvalues, and X 2 - It may be expected that approximations to the 
derivatives of A with respect to and X 2 » which are available from eqs. (14) at the convergence 
of the iterative scheme of eq. (8), can be used in eqs. (31) and (32). However, this is not advisable 
because these approximations do not necessarily converge to exact values even though the cor- 
rections in the iterative scheme converge to Newton corrections. This behavior is demonstrated in 
Table 3 for typical elements of the system matrix of example 1. Table 3 demonstrates that the de- 
rivatives of A at the convergence of the iterative scheme of eq. (9) deviate significantly from exact 
values, unless the initial guess is extremely good. Large deviations were observed in the case of 
example 2 also. As pointed out by Dennis and More (1977), this behavior is not unexpected for 
certain quasi- Newton methods including Broyden's method. Hence, it is recommended that the 
derivatives of A be re-evaluated at the solution for use in sensitivity analysis. 


Concluding Remarks 

In this paper, the complex transcendental eigenproblem with pairs of real eigenvalues is considered. 
An improved computational method for obtaining the eigenvalues of this eigenproblem is pre- 
sented. The method is based on determinant iteration and achieves a high rate of convergence that 
is close to that of Newton's method with the same cost per iteration as the generalized secant and 
Broyden's methods. The method utilizes the philosophy of Broyden's updates, except that the 
updates are applied to the system matrix derivatives rather than the Jacobian matrix as in Broyden's 
method. Computational algorithms are also presented for calculating the sensitivities of these 
eigenvalues and the corresponding eigenvectors with respect to design parameters. One algorithm 
computes the sensitivities of the eigenvalues only, while the second algorithm computes the sensi- 
tivities of the eigenvalues as well as the eigenvectors in an integrated manner. The algorithms pre- 
sented are verified by applying them to example problems. These algorithms are expected to prove 
useful in aeroelastic sensitivity analysis and optimization procedures. 
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Table 1. - COMPARISON OF RATE OF CONVERGENCE - Example 1. 
( Eigenvalue-pair of interest is 2] = 0.5 and ?. 2 = 1.5 ) 


Initial 

guess 

Number of Iterations for Convergence 


^2 

Newton's 

Present 

Broy den's 

method 

method 

method 

0.4 

1.5 

3 

4 

4 

0.7 

1.5 

3 

4 

6 

0.3 

1.3 

3 

5 

6 

0.3 

1.7 

3 

5 

6 

0.7 

1.3 

3 

5 

6 

0.7 

1.7 

3 

5 

6 

0.9 

1.9 

4 

6 

10 

0.9 

1.1 

4 

6 

10 

0.4 

1.9 

4 

5 

13 

0.4 

1.1 

4 

5 

13 
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Table 2. - COMPARISON OF RATE OF CONVERGENCE - Example 2. 
( Eigenvalue-pair of interest is = 0.5 and / 2 = 0.1 ) 


Initial guess 

Number of Iterations for Convergence 

^1 

k 2 

Newton's 

method 

Present 

method 

Broyden's 

method 

0.51 

0.11 

5 

5 

19 

0.54 

0.14 

9 

10 

No Conver- 

0.60 

0.20 

14 

20 

gence 

No Conver- 

0.58 

0.09 

6 

6 

gence 

21 

0.49 

0.11 

5 

5 

14 

0.49 

0.12 

6 

6 

38 

0.48 

0.11 

4 

4 

8 


Table 3. - ACCURACY OF MATRIX DERIVATIVES AT CONVERGENCE - Example 1. 
( Eigenvalue-pair of interest is 2] = 0.5 and ./ 2 = 1.5 ) 


^ 11^1 -0 ( -2 

a;., j exact a/j j exact ” 


Initia] 

guess 

Derivative at Convergence 

'l 

'•2 

da ]2 ldk j 

5^2/ • 1 

0.4 

1.5 

0.000 

2.000 

0.7 

1.5 

0.000 

2.000 

0.3 

1.3 

0.104 

1.947 

0.3 

1.7 

-0.104 

1.947 

0.7 

1.3 

-0.104 

1.947 

0.7 

1.7 

0.104 

1.947 

0.9 

1.9 

0.404 

1.797 

0.9 

1.1 

-0.404 

1.797 

0.4 

1.9 

-0.085 

1.961 

0.4 

1.1 

0.085 

1.961 
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